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ABSTRACT 



Context. The high-frequency quasi-periodic oscillations (HF QPOs) in neutron star and stellar-mass black hole X-ray 
binaries may be the result of a resonance between the radial and vertical epicyclic oscillations in strong gravity. 
Aims. In this paper we investigate the resonant coupling between the epicyclic modes in a torus in a strong gravitational 
field. 

Methods. We perform numerical simulations of axisymmetric constant angular momentum tori in the pseudo-Newtonian 
potential. The epicyclic motion is excited by adding a constant radial velocity to the torus. 

Results. We verify that slender tori perform epicyclic motions at the frequencies of free particles, but the epicyclic 
frequencies decrease as the tori grow thicker. More importantly, and in contrast to previous numerical studies, we do 
not find a coupling between the radial and vertical epicyclic motions. The appearance of other modes than the radial 
epicyclic motion in our simulations is rather due to small numerical deviations from exact equilibrium in the initial 
state of our torus. 

Conclusions. We find that there is no pressure coupling between the two axisymmetric epicyclic modes as long as the 
torus is symmetric with respect to the equatorial plane. However we also find that there are other modes in the disc 
that may be more attractive for explaining the HF QPOs. 

Key words, accretion: accretion discs - black hole physics - oscillations 



1. Introduction 

Observations by Rossi-XTE have revealed that some accret- 
ing neutron stars (NS) and stellar-mass black holes (BH) 
show quasi-periodic oscillations in their X-ray fluxes with 
frequencies of 100 - 1 300 Hz (see van der Klis 2000] for 
a review). Some X-ray sources are simultaneously vary- 
ing on two frequencies that form a ratio of 3:2. For BH 
QPOs this was first noticed by Abramowicz & Kluzniak 
(|200l|) (see also Strohmaye rEOOTl McClintock & Remillard 
2004). Abramowicz et al. (2003) argued that this is also 
the case for NS QPOs: although the ratio is varying, it 
clearly distinguishes the same 3:2 value (see Belloni et al. 
2005 for a criticism). The mechanism that is responsible 
for this phenomenon is not fully understood. It has been 
pointed out by Kluzniak & Abramowicz (|2000|) that the 
rational frequency ratio could be formed naturally in a res- 
onance between two modes of accretion disc oscillations. 
Radial and vertical epicyclic oscillations in an accretion 
disc have been the most often discussed possibility (e.g. 
Abramowicz & Kluzniak 2004, or a collection of review ar- 
ticles in Abramowicz 2005) . In Newtonian gravity with the 
1 jr potential the epicyclic frequencies lur and lu z are both 
equal to the Keplerian frequency £!k- This is no longer the 
case in the strong gravitational field close to a compact ob- 
ject, where lu z > ujr, which allows the two frequencies to 
assume a 3:2 commensurability at some location. 

Several authors have numerically studied oscillations of 
thick accretion discs around black holes and neutron stars 
in the QPO context in both the Newtonian (Rubio-Herrera 
& Lee I2005a .b) and the relativistic regime (Rezzolla et al. 



2003a,b, Montero et al. 120031 Za notti et al. l2005jl . Recently, 
Lee, Abramowicz & Kluzniak ( |2004[) have performed nu- 
merical simulations of pseudo-Newtonian slender tori in or- 
der to study the response of the torus to impulsive and 
periodic perturbations associated with the spin of the neu- 
tron star. They reported that a vertical epicyclic motion 
was excited in the torus due to the non-linear coupling be- 
tween the epicyclic modes when the applied perturbation 
was purely radial, and found its response to a periodic per- 
turbation to be dependent on the difference between the 
spin frequency and the frequency difference between the two 
epicyclic modes (see also Kluzniak et al. 12004 and Kluzniak, 
Abramowicz & Lee l2004p . 

In this paper we present a numerical study of oscilla- 
tory modes of slender tori in equilibrium that are given a 
small impulsive radial perturbation. We start in Sect. 2 by 
writing down the hydrodynamic equations, from which we 
construct our equilibrium tori, and then we describe the 
numerical method that we use to follow the time evolution 
of the tori. The results of our numerical simulations are de- 
scribed in Sect. 3, while Sect. 4 is devoted to a discussion 
of the results and our conclusions. 



2. Model 

2.1. Hydrodynamics 

We use the equations of nonrelativistic ideal hydrodynamics 
to describe the oscillations of a torus around a black hole. 
Ignoring the effects of radiative transport and self-gravity, 
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the equations take the form 



dp 



V • (pv) = 0, 



dt 
dv 

p— + (pv • V)v 



dt 
d(pe) 
dt 



-VP-pV$, 
V • (pev) = -PV • v, 



(1) 
(2) 
(3) 



where p and e are the density of mass and the specific in- 
ternal energy, v is the fluid velocity, <I> the gravitational 
potential and P is the pressure, which is given by the equa- 
tion of state for an ideal gas, P = pe(j — 1), with 7 = 4/3. 

The most important general relativistic effects in the 
Schwarzschild metric are mimicked by using the pseudo- 
Newtonian potential (Paczyhski & Wiita [TMD)) 



$ = - 



GM 



(r 



(4) 



where r is the spherical radius and r g = 2GM/c 2 is 
the Schwarzschild radius of the black hole. In this poten- 
tial, the epicyclic frequencies have qualitatively the same 
behaviour as in the Schwarzschild spacetime, i.e. ujr < 
u z = f^K (see Fig. [T]), and are explicitly given by wj; = 
2nv R = y/GM (R - 3r g )/ v / R(R - r g ) 3 and iv z = 2ttv z = 
y/GM/R/(R — r g ). In the Schwarzschild metric, the fre- 
quencies assume the 3:2 ratio at R = 10.8 GM/c 2 , while in 
the Paczyhski & Wiita potential this resonance occurs at 
R = 9.2 GM/c 2 . 

We construct a non-accreting torus in hydrostatic equi- 
librium with constant specific angular momentum I = Iq as 
our initial state. For convenience we will adopt geometric 
units, where G = c = 1, for the rest of the paper. The only 
non-zero velocity component is then the azimuthal compo- 
nent, which in cylindrical coordinates {P, 4>, z} is given by 

— CIR. The equilibrium configuration is then governed 
by the time-independent Euler equation 



-VP = v$ - n 2 R. 
p 



(5) 



For a barotropic fluid, in which Q — £l(R), we can define 
a rotational potential <E> ro t = / fl 2 RdR = / l 2 /R 3 dR. In 
addition we use a polytropic equation of state, P = Kp 1 , 
where K denotes the polytropic constant. The Euler equa- 
tion can then be integrated to give an expression for the 
equipotentials W that determine the shape of the torus: 
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P 



(7 - 1) P 



off 



C = W = const. 



(G) 



Here $ e ff = <& + <J> rot is the effective potential and C an inte- 
gration constant defining the surface of the torus. For par- 
ticular values of the specific angular momentum, an equilib- 
rium torus can now be constructed by filling the appropri- 
ate equipotential surfaces with a rotating fluid. An example 
of the meridional cross-section of a slender (in the sense that 
its radial extent is much smaller than its distance from the 
central object) torus is shown in Fig. [2] 

2.2. Numerical simulations 

The numerical simulations in this paper are performed us- 
ing the two-dimensional ZEUS-2D code by Stone & Norman 
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Fig. 1. Radial (ojr) and vertical (ui z — f2x) epicyclic fre- 
quencies in the Paczyhski & Wiita potential, and the cor- 
responding radial (lower dashed curve) and vertical (up- 
per dashed curve) epicyclic frequencies in the Schwarzschild 
metric, for a M = 10 M Q black hole. 



0.10 



0.05 



N 



0.00 



-0.05 




-0.10 h ,,,,,,,,,,,,,,,,,,,,,,, - 

7.90 7.95 8.00 8.05 8.10 
R/M 

Fig. 2. Density contours in a meridional cross-section of a 
slender torus at R — 8 M with radial extent of approxi- 
mately 0.15 M. 



(fT992|) that solves Eqs. JTJ-©. ZEUS is a time-explicit, 
Eulerian finite difference code that uses a staggered grid. 
Shock waves are handled by an artificial viscosity tensor. 
The magnetic field, self-gravity and radiation transport in 
ZEUS were switched off during our simulations. 

We study two different groups of models (Tab. Q] ). 
Models A1-A4 represent slender tori with a radial extent 
of approximately 0.15 M at different radial positions (see 
Fig. [2]). Models B1-B4 make up a series of thicker (non- 
slender) tori. We perturb the tori by adding a velocity field 
that is constant in space at t = 0. The perturbation is 
always subsonic, corresponding to a Mach number at the 
torus center of mass of 0.3. 

We use a uniform grid in cylindrical coordinates in all 
our simulations and the resolutions of the different models 
are specified in Tab. Q] 
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Fig. 3. Radial (top) and vertical (bottom) oscillations of a 
torus centered at R — 9.2 M (Model A3) where the two 
epicyclic modes are in a 3:2 ratio. The time is given in 
units of the Keplerian period Tk defined at the centre of 
the torus. The torus performs a radial oscillation with an 
amplitude of 0.01 M . The vertical oscillation excited here 
is seven orders of magnitude smaller than the radial oscil- 
lation. 
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Fig. 4. Same as the bottom panel in Fig. [3l but using a 
finer grid (Model A4). Note the decrease of the amplitude 
of the oscillation compared to Fig. [3l 



3. Results 

For the data analysis we look at some of the characteris- 
tic quantities as functions of time, such as the total kinetic 
energy, the mean of the square of the density (as was pre- 
viously done by e.g. Zanotti et al. 120051 and Blaes et al. 
2006a) , and (following Lee et al. I2004p the radial and ver- 
tical positions of the centre of mass of the torus, as well as 
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Fig. 5. Top: Power spectrum of (p 2 ) of a radially perturbed 
slender torus (Model A2). Frequencies are given in units 
of f^K, the Keplerian angular velocity at the centre of the 
torus. The left peak at approximately 0.58 J1k corresponds 
to the radial epicyclic motion (see Tab.QJ. The other main 
peaks are located at about 0.86 11k and 1.5 £Ik- The vertical 
lines denote the frequencies of the radial epicyclic mode (a), 
the plus- mode (b) and the breathing mode (c). Bottom: 
Same as the top panel, but for an unperturbed torus. The 
peak belonging to the radial motion is absent, but the other 
two peaks remain with approximately the same power, as 
in the top panel. 



snapshots of the density, internal energy and velocity field 
inside the torus. 



3.1. Slender torus oscillations 

We consider a slender torus (Models A1-A4) and perturb 
it slightly by adding a purely radial constant velocity field 
at t = 0. The torus then performs radial oscillations at 
the radial epicyclic frequency of a test particle in a circular 
orbit at the centre of the torus (Fig. [3]). We also find a small 
oscillation in the vertical direction at the vertical epicyclic 
frequency of a test particle (see Tab. Q] for the frequencies 
measured from the simulations), but the amplitude of this 
motion is more than six orders of magnitude smaller than 
that of the radial epicyclic motion. The small amplitude of 
the vertical oscillation is independent of both the amplitude 
of the radial epicyclic motion and the radial position of the 
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Fig. 7. Snapshots of the pv inside the torus (Model A2) after substraction of the radial epicyclic motion. The snapshots 
cover a period of the breathing mode, T — (1/1.54) Tk, where Tk again denotes the Keplerian period at the centre of 
the torus. The residual pw resembles that of the breathing mode (see Fig. 1 in Blaes et al. 12006b ). 



centre of the torus. Moreover, it decreases with increasing 
numerical resolution, as shown in Fig. [4] Our conclusion is 
therefore that it is not the result of a resonance between 
the two epicyclic modes. Rather it is excited by numerical 
noise in the simulations. 

The power spectra of the mean of the square of the 
density and the total vertical kinetic energy for Model A2 
can be seen in Figs. [5] and [SJ The (p 2 ) power spectrum is 
not sensitive to incompressible modes, such as the verti- 
cal epicyclic mode, but the radial epicyclic mode is weakly 
compressible because of the cylindrical geometry, and thus 
we can see it in Fig. [S] It corresponds to the left peak at 
ojr = 0.57 S!k- Beyond that, there are two other peaks 
in the power spectrum of (p 2 ). These peaks appear inde- 
pendently of whether the torus has been perturbed or not 
(compare the top and bottom panels of Fig. [5]). The strong 
peak at 1.5 fix is present in the spectrum at this frequency 
regardless of the location of the torus, while the less promi- 
nent peak at approximately 0.86 Ox slightly changes its 
position depending on the location of the torus. 

Blaes et al. (|2006r ) have calculated the lowest-order 
slender torus modes. Two of these modes, the incompress- 



ible "plus- mode" and the acoustic "breathing" mode have 
frequencies of 0.87 Ox and 1.54 Ox for the torus of Model 
A2 (see their Tab. 2 for explicit expressions). These fre- 
quencies match well with the peaks at 0.86 Ox and 1.5 J7k 
in Fig. G3 

The 1.5 £1k peak also appears in the power spectrum of 
the total vertical kinetic energy, as seen in Fig. [S] We also 
find four minor peaks in the power spectrum of the vertical 
kinetic energy at 0.45 17k, 0.75 fix, 1-05 fix, and 1.2 Ok- 
The peaks at 0.45 fix and 1.05 Ok are absent in the case 
of an unperturbed torus (bottom panel of Fig. [6]). None of 
the minor peaks have any corresponding peaks in the power 
spectrum of (p 2 ), and they cannot be identified with any of 
the slender torus modes that Blaes et al. (2006a) calculated. 
An intriguing feature of Fig. [S] is the lack of power at the 
frequency of the radial epicyclic mode, which suggests that 
this mode has a negative influence on the vertical motion. 

Fig. [7] shows four snapshots of pv inside the torus 
(Model A2) after substraction of the mean radial motion 
of the torus. They cover one period of the breathing mode 
2-7r/(1.54 fix)- The velocity pattern is dominated by the 
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Table 1. Properties of the models. From left to right the columns give the name of the model, grid resolution, radial 
(i?-grid) and vertical (z-grid) boundaries of the grid, the inner (i?i n ), central (R c ), and outer (i? ut) radii, the initial 
velocity kick in terms of the Mach number M at the centre of the torus, and radial lur and vertical lo z epicyclic frequencies 
of the torus centre in units of fix- 
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Fig. 6. Top: Power spectrum of the total vertical kinetic 
energy of a radially perturbed slender torus centered at 
R = 8 M (Model A2). Frequencies are given in units of 
Qk- Note the major peak at 1.5 Qk and the smaller peaks 
at 0.45 Cl K , 0.75 K , 1.05 fi K , and 1.2 tt K . The vertical 
lines denote the frequencies of the vertical epicyclic mode 
(a) and the breathing mode (b) calculated for this Model. 
Bottom: Same as the top panel, but for an unperturbed 
torus. The 0.45 Ok and 1.05 f^K peaks are missing in this 
case 



breathing mode, though we do see radial flows that match 
the plus-mode too. 
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Fig. 8. Radial ur and vertical uj z epicyclic frequencies of 
the torus as a function of its width assuming that the torus 
is centered at 9.2 M. 



3.2. Epicyclic oscillations in thick tori 

In order to excite both the epicyclic modes, we applied 
a perturbation in both the radial and the vertical direc- 
tions. The resulting epicyclic frequencies for four different 
tori (Models B1-B4) are presented in Tab. [T] For an in- 
finitely slender torus, the frequencies are consistent with 
the epicyclic frequencies of a test particle at the torus cen- 
tre, lor — 0.666 f^K and u> z — £!k- The frequencies decrease 
as the torus grows thicker, as seen from TabQ] and Fig. 
[5] A similar trend was reported by Rubio-Herrera & Lee 
(|2005b l. and Blaes et al. (|2006b ). where it has been ex- 
plored analytically for slightly non-slender Newtonian tori. 
It is worth mentioning that even when we excite both the 
epicyclic modes, we do not see any exchange of energy be- 
tween them (see Fig. [9]). 



4. Discussion and Conclusions 

4.1. Epicyclic oscillations 

We have performed two-dimensional numerical simulations 
of oscillations of constant angular momentum tori in the 
Paczyhski & Wiita potential. The tori were given an impul- 
sive radial perturbation and were allowed to evolve in time 
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Fig. 9. Radial (top) and vertical (bottom) oscillations of 
a torus centered at R = 9.2 M (Model A3) when both the 
modes are excited. Note that there is no exchange of energy 
between the two modes. 



afterwards. The tori are then describing radial epicyclic mo- 
tions. For a sufficiently slender torus the frequency of this 
oscillation agrees with that of a test particle. As the torus 
grows thicker we find that the epicyclic frequencies signifi- 
cantly decrease. This is in accord with previous numerical 
results by Rubio-Herrera & Lee p005b ). as well as recent 
analytic work where the expressions for the (negative) cor- 
rections to the epicyclic frequencies were derived using a 
perturbative method accurate to second-order in the ex- 
tent of the torus (Blaes et al. 12006b ). 

One of the main objectives of this work is to investi- 
gate the modes that may be excited by a radial epicyclic 
oscillation in the torus. In particular, we look for the reso- 
nant coupling between the radial and the vertical epicyclic 
modes, which Kluzniak & Abramowicz suggested plays a 
significant role in the high-frequency QPOs. We find that 
the centre of mass of the torus performs a vertical oscilla- 
tion of a very low amplitude, but it is independent of the 
amplitude of the radial epicyclic oscillation and the loca- 
tion of the torus, unlike what one would expect if it was due 
to a resonance between the two epicyclic modes. Rather it 
is likely that the vertical epicyclic oscillation is excited by 
numerical noise in the simulations. A coupling between the 
radial and vertical epicyclic modes in the 3:2 parametric 
resonance may only happen in a limited range of the pa- 
rameter space, the Arnold tongue (see e.g. Arnold 119781 
Arnold 1983). This range is not known a priori, but one 
relevant parameter is the initial vertical amplitude. 

Lee et al. (|2004p found that a small vertical epicyclic 
motion was excited by the radial epicyclic motion in their 
simulations of an axisymmetric torus. We believe that the 



reason for this is that their initial torus was not perfectly 
symmetric around the z = plane, and thus the initial 
amplitude of the vertical oscillation was not zero. We did 
not find that any vertical oscillations were excited in our 
simulations, in which the initial amplitude of the vertical 
motion was zero - thus in our case there was not any cou- 
pling between the radial and vertical epicyclic modes. This 
is hardly surprising: conservation of momentum prohibits 
the excitation of a vertical epicyclic motion of the centre 
of the torus, if it is initially at rest in the equatorial plane, 
and no external vertical forces are at play. 

In order to find the Arnold tongue empirically one 
should perform many simulations similar to these described 
in our paper, but with different initial vertical amplitudes, 
and with different ratios (close to 3:2) of the radial and ver- 
tical epicyclic frequencies. This will give the range in the 
parameter space in which the mode coupling for axisym- 
metric oscillations may occur. 

A more promising approach might be to consider non- 
axisymmetric oscillations. The simplest example of nonax- 
isymmetric vertical motion is an m = 1 precessing warp in 
an accretion disc. Such a warp will resonantly drive a radial 
epicyclic motion (Papaloizou & Pringle[l983, Papaloizou & 
Lin [1995). A warp can be produced, for instance, by the 
radiation pressure from a central radiation source (Pringlc 
1996). The role of m = 1 warps for the HF QPOs has been 
discussed by Kato (PDMI2Uo1))) . 

Another possibility is that the two modes can be cou- 
pled by the local fluctuations that are generated by the 
turbulence in the disc (note that the turbulence is non- 
axisymmetric). It is today widely accepted that accretion 
discs become turbulent because of the magnetorotational 
instability, that was introduced by Balbus & Hawley (1991). 
This turbulence generates the radial transport of angular 
momentum that drives the accretion and heats the disc, but 
its role in other aspects of the dynamics of accretion discs 
has hardly been explored though. Vio et al. (|2006|) have 
recently studied how stochastic processes can produce res- 
onances between the epicyclic modes in an almost Keplerian 
disc. 



4.2. Other oscillatory modes 

A problem with attributing the HF QPOs to a resonance 
between the two epicyclic modes is that they only occur 
at a special resonance radius, and thus require a certain 
amount of fine-tuning. However our tori contain more os- 
cillatory hydrodynamical modes, and the power spectra of 
(p 2 ) reveal two significant peaks, both of whom occur re- 
gardless of whether there is an initial perturbatiorQ. Firstly, 
there is a peak at 1.5 S7k, which is also displayed in the 
power spectrum of the vertical energy. This peak appears 
at approximately the same frequency in all the simulations, 
no matter what is the location of the torus. It is of interest 
since its frequency forms a ratio of 3:2 with the Keplerian 
frequency 17k, which is also the vertical epicyclic frequency. 
It can perhaps be identified as the lowest-order acoustic 
wave mode (see Blaes et al. 12006b ). which has a frequency 
close to 1.5 Ok, but with some dependence on the radial lo- 
cation of the torus. It was suggested by Blaes et al. (p006a) 



1 These modes are excited in our simulations because the equi- 
librium torus of Sect. 2 is not any longer in perfect equilibrium 
when it is implemented on our grid. 
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that this mode and the vertical epicyclic mode might be 
responsible for the 3:2 commensurability observed in the 
HF QPOs. 

The other peak at 0.86 Qk agrees well with the fre- 
quency of a surface gravity wave mode that Blaes and 
his collaborators labeled as the plus-mode. This oscillatory 
mode was found to be excited in previous numerical stud- 
ies with either localized (Rubio-Herrera & Lee |2005h,,b), 
or global (Rez zolla e t al. [2003k .b. Montero et al. 120041 and 
Zanotti et al. 1 200 5]) external perturbations, where it was 
attributed to one of the set of acoustic p-modes of oscil- 
lations. Such identification was based on numerical calcu- 
lations of eigenfunctions and eigenfrequencies of vertically 
integrated tori, carried out by Rezzolla et al. 2003b, and it 
was clarified in Blaes et al. (12006a ) (see their discussions 
in Sects. 4.4 and 5.1) that this in general incompressible 
mode may indeed be related to an acoustic p-mode in the 
two-dimensional case of a vertically integrated torus. The 
frequencies of the p-modes were found to appear in an ap- 
proximate harmonic sequence 2 : 3 : 4..., in which the first 
frequency is consistent with the radial epicyclic frequency 
of the torus, the second is the plus-mode frequency, and 
others correspond to the higher-order p-mode frequencies. 
The plus-mode and the radial epicyclic mode that appear 
in a 3:2 ratio were proposed by Rezzolla et al. (2003a) as 
an explanation of the black hole HF QPOs. 
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